###############################################
## BAR PLOT FUNCTIONS
###############################################

## bar plots for sequences with zero reads
make_zero_read_qc_reads_plot <- function(df){
## count rows 
  plot_df <- df %>%
  select(lab, target_site, candidate_off_target_site) %>% 
  count(lab, target_site)
  
  ## plot 
  ggplotly(ggplot(plot_df) + 
    geom_bar(aes(x = target_site, y = n, fill = lab), 
             position = position_dodge(preserve = "single"), stat = "identity") +
    scale_y_continuous(labels = comma) +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    labs(x = "Target Site", y = "Number of Sequences with 0 reads", 
         fill = "Lab") + scale_fill_manual(values=cbPalette))
}  

make_read_cutoff_bar_plot <- function(df, read_count_value){
## count rows 
  plot_df <- df %>%
  filter(post_selection_read_count >= read_count_value) %>% 
  select(lab, target_site, candidate_off_target_site) %>% 
  count(lab, target_site)
  
  plot_title <- paste("Number of Barcodes with at least", {{read_count_value}}, "Reads")
  
  ## plot 
  ggplotly(ggplot(plot_df) + 
    geom_bar(aes(x = target_site, y = n, fill = lab), 
             position = position_dodge(preserve = "single"), stat = "identity") +
    scale_y_continuous(labels = comma) +
    theme_bw() + 
    theme(axis.text.x = element_text(angle = 90)) + 
    labs(x = "Target Site", y = "Count", 
         fill = "Lab", title = plot_title) + scale_fill_manual(values=cbPalette))
}  

make_outlier_count_plots <- function(df){
  
  plot_df <- df %>% 
    select(post_selection_read_count, candidate_off_target_site, lab, target_site,
           CI_range) %>% 
    count(target_site, lab, CI_range) %>%
    rename(num_barcodes = n)
  
  ggplotly(ggplot(plot_df) +
      geom_bar(aes(x = target_site, 
                   y = num_barcodes,
                   fill = CI_range),
               position = position_dodge(preserve = "single"), 
               stat = "identity"
               ) +
          scale_y_continuous(labels = comma) +
      ## Split plot by summary metric
      facet_wrap(~lab, scales = "free_y") +
      theme_bw() +
      theme(axis.text.x = element_text(hjust = 0)) +
    labs(x = "Target Site", y = "Number of Barcodes", fill = "95% CI", 
         title = "Number of Barcodes Above and Below 95% Confidence Interval") +
      scale_fill_manual(values=cbPalette)) %>% 
    layout(legend = list(orientation = "h", x = 0.4, y = -0.2))
}

make_outlier_read_count_plot <- function(df, targ_site){
  targ_site_plot_df <- df %>% 
    filter(target_site == {{targ_site}})
  
  ggplotly(ggplot(targ_site_plot_df) +
    geom_bar(aes(x = post_selection_read_count, 
                 y = n,
                 fill = CI_range),
             position = position_dodge(preserve = "single"), 
             stat = "identity"
             ) +
        scale_y_continuous(labels = comma) +
    ## Split plot by summary metric
    facet_wrap(~lab, scales = "free_y") +
    theme_bw() +
    theme(axis.text.x = element_text(hjust = 0)) +
  labs(x = "Read Count", y = "Count", fill = "95% CI", 
       title = "Number of Barcodes") +
    scale_fill_manual(values=cbPalette)) %>% 
  layout(legend = list(orientation = "h", x = 0.4, y = -0.2))
}

make_outlier_read_count_table <- function(df, targ_site){
  (targ_site_plot_df <- df %>% 
    filter(target_site == {{targ_site}}) %>% 
    pivot_wider(names_from = "CI_range", values_from = "n", values_fill = 0) %>% 
    rename(Number.Barcodes.Above.95CI = above, 
           Number.Barcodes.Below.95CI = below,
           Target.Site = target_site,
           Lab = lab, 
           Read.Count = post_selection_read_count) %>% 
     arrange(Read.Count) %>%
    kbl() %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")) %>% 
    kable_styling(font_size = 10))
}

########################################################################
#  READ COUNT SUMMARIES
########################################################################
## get list of files in results directory
QC_reads_file_list <- as.list(list.files(path = here("data/Vikram_QC/report_input"), 
                           pattern = "QC.csv",
                           include.dirs = TRUE, full.names = TRUE, 
                           recursive = TRUE))

## strip the path portion of the file names
QC_reads_file_names <- str_remove(QC_reads_file_list, 
                                  paste0(here("data/Vikram_QC/report_input"),"/"))

## organize the file list to read them into a df
QC_reads_lst <- set_names(QC_reads_file_list, QC_reads_file_names)

## make lists into a data frame
QC_reads_metadata_df <- tibble(QC_reads_file = unlist(QC_reads_file_names)) %>% 
  separate(QC_reads_file, c("lab","target_site"), 
           sep = "_", remove = FALSE) %>%
  mutate(target_site = str_remove(target_site, "QC.tsv"))

## read in files & annotate with the metadata
QC_reads_df <- QC_reads_lst %>%
  map_dfr(read_csv, .id = "QC_reads_file") %>% 
  left_join(QC_reads_metadata_df) 

## create a df that holds lib size info
target_site <- c('EMX1', 'FANCF', 'RNF2') 
lib_size <- c(51743, 29207, 26102) ## num rows of lib.txt files
lib_size_df <- data.frame(target_site, lib_size)

QC_reads_df <- QC_reads_df %>% 
  left_join(lib_size_df)

########################################################################
# SUMMARY OUTPUT
########################################################################

## get list of files in results directory
qc_file_list <- as.list(list.files(path = here("data/Vikram_QC/report_input"),
                           pattern = "QC.tsv",
                           include.dirs = TRUE, full.names = TRUE,
                           recursive = TRUE))

## strip the path portion of the file names
qc_file_names <- str_remove(qc_file_list, paste0(here("data/Vikram_QC/report_input"),"/"))

## organize the file list to read them into a df
qc_df_lst <- set_names(qc_file_list, qc_file_names)

## make lists into a data frame
qc_metadata_df <- tibble(qc_file = unlist(qc_file_names)) %>%
  separate(qc_file, c("lab","target_site"),
           sep = "_", remove = FALSE) %>%
  mutate(target_site = str_remove(target_site, "QC.sv"))

## read in files & annotate with the metadata
qc_df <- qc_df_lst %>%
  map_df(read_tsv, col_names = c("var","number_of_seqs"), .id = "sample") %>%
  mutate(var= str_remove(var, "number_of_")) %>%
  separate(sample,c("lab","target_site"), sep = "_", remove = FALSE) %>%
  select(-sample) %>%
  mutate(target_site = str_remove(target_site, ".QC.sv"))

########################################################################
# COLOR PALETTE 
########################################################################
## establish color palette for plots
cbPalette <- c("#000000", 
               "#E69F00", 
               "#56B4E9", 
               "#009E73", 
               "#F0E442", 
               "#0072B2", 
               "#D55E00", 
               "#CC79A7", 
               "#999999", 
               "#20F03B",
               "#FEC44F", 
               "#D95F0E", 
               "#756BB1", 
               "#FF8A33", 
               "#D7BE07", 
               "#A807D7")

########################################################################
# CLEAR OBJECTS 
########################################################################
rm(QC_reads_file_list)
rm(QC_reads_file_names)
rm(QC_reads_lst)
rm(qc_file_list)
rm(qc_file_names)
rm(qc_df_lst)
rm(QC_reads_metadata_df)
rm(qc_metadata_df)

## df of all zero read sequences
zero_count_df <- QC_reads_df %>%
  filter(post_selection_read_count == 0)

## plot how many zero-read count seqs per target site and lab
make_zero_read_qc_reads_plot (zero_count_df)

###############################################################
## sanity check - calculate % library members not present
## against Vikram's QC script output
  
## create sanity check df and perform calculations
zero_count_sanity_check_df <- zero_count_df %>% 
  select(lab, target_site, candidate_off_target_site) %>% 
  count(lab, target_site) %>% 
  left_join(lib_size_df) %>%
  group_by(target_site, lab) %>% 
  mutate(pct_lib_mems_not_present = round((n/lib_size)*100, digits = 3)) %>%
  rename(num_seqs_with_zero_reads = n) 

## compare these values to the values calculated in Virkam's QC table
(zero_count_qc_check <- qc_df %>% 
  filter(var == "pct_lib_members_not_present") %>% 
  pivot_wider(names_from = var, values_from = number_of_seqs, values_fill = NA) %>%
  right_join(zero_count_sanity_check_df) %>% 
    arrange(target_site) %>%
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")) %>% 
    kable_styling(font_size = 10))

## plot 
make_read_cutoff_bar_plot(QC_reads_df, 6)

## table
(cutoff_count_df <- QC_reads_df %>%
  filter(post_selection_read_count >= 6) %>% 
  select(lab, target_site, candidate_off_target_site) %>% 
  count(lab, target_site) %>% 
  rename(num_seqs_w_at_least_6_reads = n) %>% 
    arrange(target_site) %>%
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")))

## create new df with all of these calculations
stats_df <- QC_reads_df %>% 
  group_by(target_site, lab) %>% 
  mutate(average_read_count = mean(post_selection_read_count), 
         std_dev_read_count = sd(post_selection_read_count),
         two_sigma_minus = average_read_count-(2*std_dev_read_count),
         two_sigma_plus = average_read_count+(2*std_dev_read_count),
         num_barcodes_less_95CI = sum(post_selection_read_count < two_sigma_minus),
         num_barcodes_more_95CI = sum(post_selection_read_count > two_sigma_plus),
         pct_barcodes_less_95CI = paste0(round((num_barcodes_less_95CI/lib_size)*100, digits = 1), "%"),
         pct_barcodes_more_95CI = paste0(round((num_barcodes_more_95CI/lib_size)*100, digits = 1), "%")) %>% 
  mutate(CI_range = case_when(post_selection_read_count < two_sigma_minus ~ "below",
                              post_selection_read_count > two_sigma_plus ~ "above",
                              post_selection_read_count > two_sigma_minus & post_selection_read_count < two_sigma_plus ~ "within"))
                              

## downsize df and show summary values
(distinct(stats_df %>% 
    select(lab, target_site, lib_size, average_read_count, std_dev_read_count, 
           two_sigma_minus, two_sigma_plus, num_barcodes_less_95CI, num_barcodes_more_95CI, 
           pct_barcodes_less_95CI, pct_barcodes_more_95CI)) %>% 
    arrange(target_site) %>% 
  kbl() %>%
  kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")) %>% 
    kable_styling(font_size = 10))

## reduce df to necessary columns, pivot to plot
outliers_df <- stats_df %>% 
  filter(!CI_range == "within")
  
make_outlier_count_plots(outliers_df)

(outlier_plot_df <- outliers_df %>% 
    select(post_selection_read_count, candidate_off_target_site, lab, target_site,
           CI_range) %>% 
    count(target_site, lab, CI_range) %>%
    rename(num_barcodes = n) %>% 
    pivot_wider(names_from = "CI_range", values_from = "num_barcodes", values_fill = NA) %>% 
    rename(Number.Barcodes.Above.95CI = above, 
           Number.Barcodes.Below.95CI = below,
           Target.Site = target_site,
           Lab = lab) %>% 
    kbl() %>%
    kable_styling(bootstrap_options = c("striped", "hover", "condensed", "responsive", html_font = "Cambria")))